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NUMERICAL METHODS FOR METEOROLOGY AND CLIMATOLOGY 


E. Isaacson, D. Marchesin, G. Zwas, Courant Inst of Math' I Sciences , 
New York University 


ABSTRACT 

Efficient numerical methods for long term 
weather forecasting are developed. One 
implicit and one explicit scheme are compared 
as to accuracy. 

INTRODUCTION 


It has been estimated that about 30% of the error 
in meteorological forecasts is attributable to the in- 
accuracy of numerical schemes currently used. In 
order to minimize this error, we have developed and 
tested two new numerical methods. These schemes are 
particularly designed to perform long term inte- 
grations for climatological purposes, since they 
conserve the global quantities: total mass, total 
energy, and total potential enstrophy of the flow. 

One of the schemes is implicit, while the other one is 
explicit. The greater inherent efficiency of the 
implicit scheme, for fine spatial grid calculations 
arises from the fact that its time step At may be an 
integer multiple of AX/U , and is not required to be 
less than AX/(C+U) as is the case for explicit 
schemes, where AX is the smallest spatial interval 
size, U is the maximum speed of propagation of the 
flow, while C is the much larger maximum speed of 
sound waves in the flow. In other words, a much 
larger time step is permitted for implicit schemes. 
But, in order to derive this benefit it is essential 
that the amount of computer time used to solve the 
implicit equations at any time step, be a small 
multiple of the time needed to make one time step with 
an explicit scheme. Our application of these schemes 
indicates that indeed the implicit scheme requires 
less computer time than does the explicit scheme, for 
fine grid sizes. 

In section 2, both the implicit and explicit 
methods are used to solve the shallow water equations 


183 


for a single- layered atmosphere over a smooth earth 
and comparisons of the solutions and of the computer 
time are made. In section 3 , with orography present , 
a comparison is made of the solutions produced by the 
explicit scheme both with and without the conservation 
of the global quantities. Here it is observed that 
the solution produced while conserving the global 
quantities, is smoother and more accurate. In section 
4, a brief description is given of the numerical 
schemes . 


Application of schemes 

The technique for devising an implicit scheme and 
the method for modifying both implicit and explicit 
schemes so as to conserve the global quantities were 
applied to the solution of the nonlinear shallow water 
equations for a single layered, incompressible fluid 
over a spherical earth. As initial conditions, three 
subtropical highs were centered symmetrically in each 
hemisphere on latitudes + 30°, and at longitudes 0°, 
120°E , 2 40°E . The initial wind speed is found geo- 
strophically from the initial height field. The 
initial hydrostatic "pressure" (i.e. density • height) 
rose from 1036 (corresponding to 8500 meters) to 1046 
at the center of the highs. The maximum of the 
initially geostrophic velocity is about 13 meters/sec. 
Contour levels of the initial "pressure" field are 
plotted, on a stereographic chart of the Northern 
Hemisphere, along with arrows indicating the initial 
wind speed and direction in Fig. 1. After 24 hours, 
we exhibit in Fig. 2, appropriately scaled height 
fields produced on one third of the Northern , Hemi- 
sphere with a fine spherical grid using ^60 ° in both 

longitude and latitude. Fig. 2a was produced by the 
explicit scheme, while Fig. 2b was produced by the 
implicit scheme. Note how closely the "high" cells 
resemble each other after 24 hours. The left most 
column of numbers represent longitudes, that is row 73 
corresponds to longitude 73^|P °. The rightmost 

column corresponds to the latitude nearest the equator 
and the second column from the left corresponds to the 
latitude nearest the North Pole. 
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For the explicit scheme a time step of 3 minutes 
was used, while in the implicit scheme an interval of 
15 minutes was used. The table below lists the 
computing times (CPU) required on a 60 bit machine, 
the CDC-6600, for each method with both the fine grid 
and a coarse spherical grid of 360/64 ° for latitude 
and longitude. 



explicit 

scheme 

implicit 

scheme 

revised implicit 
scheme 

coarse grid 

3 min 

9 min 

3 min 

fine grid 

24 min 

36 min 

12 min 


Computing Times on CDC-6600 

Note that the time step of 15 minutes was used in 
the implicit scheme for both the fine and coarse 
grids; whereas in the explicit scheme, 6 minutes was 
used for the coarse grid and 3 minutes was used for 
the fine grid. Both schemes produced solutions with 
comparable accuracy. Furthermore, we have observed 
that it is possible to revise the implicit scheme to 
take advantage of the special form of the hydrodynamic 
equations, so as to cut the computing time of the im- 
plicit scheme in half. This feature of the equations 
of motion is also present in the full hydrodynamic 
equations used for general circulation models. With 
this improvement the implicit scheme would run in less 
time than the explicit scheme, for fine grids. 

Orography effects 

A major source of computational errors is found 
to be the way in which current numerical methods treat 
flow over and around high and extensive mountain 
ranges, e.g. Rocky Mountains, Andes Mountains, 
Himalayan Mountains, etc. In order to eliminate this 
difficulty, Arakawa has advocated using difference 
schemes that conserve the total potential enstrophy, 
total energy, and total mass. We have found it 
possible to modify any difference scheme, so as to 
make the new scheme conserve these global quantities. 
In particular, we show the effect that this conserva- 
tive modification has when we introduced orography 
into the shallow water model previously described in 
section 2. That is, we introduced three identical 
mountain ridges centered respectively along the three 
longitudes 60°E, 180°E, and 300°E. The maximum height 
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8500 


above sea-level of the mountain ridge is 


meters , 
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at the equator. The mountain height decreases to zero 
at the poles. The initial velocity is again found 
from geostrophic balance on the sphere. The fluid 
tends to flow around the ridges, by veering to the 
pole. In Fig. 3a, we see that the solution produced 
after 24 hours without the conservation of the global 
quantities has larger velocities and steeper height 
gradients near the pole, than does the solution shown 
in Fig. 3b, which is produced with the conservation 
of total mass, total energy, and total potential 
enstrophy. In the Table below we indicate the 
relative change from the initial values of these 
global quantities, after t-24 hours and after t=48 
hours, for the explicit scheme. In the first column 
the larger relative deviations are found without the 
use of the conservative modification method; whereas 
the smaller relative deviations in the second column 
are produced with the use of the conservative modifi- 


quantity 

without 

conservative 

modification 

with 

conservative 

modification 

total mass (24 hours) 

140-10 -6 

0.2- 10" 6 


203 *10 -6 

0.5* 10~ 6 


5766 • 10 -6 

10.0-10 -6 

total mass (48 hours) 

235- 10 -6 

0.5-10 -6 


260 • 10 -6 

1.1.10 " 6 


1IWBBBB1 

16.0-10 -6 


Relative errors in total quantities. 
The Numerical Schemes 


The explicit scheme is a so-called leapfrog 
scheme in which the solution is advanced from (t-At) 
and t to the time (t+At) . The first order time 
derivatives are replaced by differences centered at 
time t while the first order spatial derivatives in 
latitude and longitude, are replaced by fourth order, 
five point, centered difference expressions. At the 
latitudes closest to the pole, the values of the 
solution determined by this difference scheme are 
smoothed by using the fast Fourier transform. This 
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program has been adapted from the ideas of Kreiss and 
Oliger [4] and Williamson and Browning [6], and is 
described in Isaacson and Stoker [5]. 

The implicit scheme is roughly speaking a Crank- 
Nicolson type scheme in which the solution is advanced 
from time t to time (t+At) . Here the fourth order 
accuracy is obtained by using a Pad*e rational fraction 
in the spatial difference operators. After clearing 
fractions, by multiplying all terms by the denominator, 
we obtain a compact three point spatial difference 
expression. The resulting simultaneous difference 
equations are simplified by writing the spatial 
operator as the product of a longitudinal factor and a 
latitudinal factor. These factors have a block tri- 
diagonal matrix representation, that is easily 
invertible. This idea and its implementation were 
proposed by A. Harten based on the methods used by 
Beam and Warming [1,2] in aerodynamics. Considerable 
effort was needed to find an efficient and stable 
factorization method. Finally, it was found necessary 
to use a fourth order Shapiro filter at each time step 
in order to maintain stability. 

The conservative modification method is described 
in general terms in Isaacson [3]. Here the work in- 
volved in modifying the solution found at time (t+At) 
is of the order of the amount of calculation used in 
one time step of the explicit scheme. 
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Fig. 2b, Height field after 24 hours - implicit scheme 
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